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ABSTRACT 

Aims. We study the accretion of dust particles of various sizes onto embedded massive gas giant planets, where we take into account 
the structure of the gas disk due to the presence of the planet. The accretion rate of solids is important for the structure of giant planets: 
it determines the growth rate of the solid core that may be present as well as their final enrichment in solids. 

Methods. We use the RODEO hydrodynamics solver to solve the flow equations for the gas, together with a particle approach for the 
dust. The solver for the particles' equations of motion is implicit with respect to the drag force, which allows us to treat the whole 
dust size spectrum. 

Results. We find that dust accretion is limited to the smallest particle sizes. The largest particles get trapped in outer mean-motion 
resonances with the planet, while particles of intermediate size are pushed away from the orbit of the planet by the density structure 
in the gas disk. Only particles smaller than approximately s mlix = 10 /im may accrete on a planet with the mass of Jupiter. For a ten 
times less massive planet s max = 100 fim. The strongly reduced accretion of dust makes it very hard to enrich a newly formed giant 
planet in solids. 

Key words. Hydrodynamics - Methods: numerical - Planets and satellites: formation 



1. Introduction 

Planets form in circumstellar disks consisting of gas and dust. 
Initially, gas and solids are well-mixed with a dust-to-gas ratio of 
1 : 100, which is similar to the interstellar value, and the dust par- 
ticles will be as small as interstellar grains (~ 1 yum). However, 
due to the pressure structure of the disk the dust particles start 
to move with respect to the gas: the vertical component of the 
stellar gravity makes the particles rain down onto the midplane 
of the disk, while the radial pressure g radient of the gas pushes 
the particles inward (Weidenschilling 1977). 

The magnitude of the velocity difference between gas and 
dust depends on the particle size: small particles couple well to 
the gas and move slower than the larger particles. These size- 
dependent velocities in the radial and vertical direction, together 
with Brownian motion and turbulence, will lead to collisions be- 
tween the dust particles of various sizes. Assuming that impact 
velocities are low the particles stick together upon colliding, af- 
ter which they continue as one large particle. This leads to a rapid 
depletion of small grains (Dullem ond & Dominikl200^) . 

However, sticking probabilities are still u nclear, especially 
for larger particles. Expe rimental results (|Blum & Muench 
fl993l lBl um & Wurm 2000) suggest that unless the relative ve- 
locities of particles are very low collisions lead to bouncing or 
even shattering of the particles. Another problem arises once 
the particles grow to sizes of approximately one meter, because 
then the time scale for drag-induced ra dial migration become s 
as short as 100 yrs for particles at 1 AU (Weidenschilling 19771) . 
These problems led to a renewed interest in planetesi mal forma- 
tion through gravita tional instability (Gold reich & Wardll 19731 
lYoudin & Shul2002HGaraud et all2004l) ~ 
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Once the planetesimals reach sizes of roughly 1 km the 
aforementioned problems cease to exist. These bodies couple not 
well enough to the gas to experience fast radial migration, and 
their gravitational influence is aiding their growth. An important 
consequence of gravity-aided growth is that the heavier bodies 
will grow fastest, which leads to a phase of oligarchic growth 
wher e a few large bodies domina te the dynamics of the system 
dKokubo & Idall998ll2^00ll200 2i This is an important epoch in 
planet formation, because it covers the phase from planetesimals 
to planets of masses comparable to Mercury or Mars in the inner 
regions of the solar system, while in the outer regio ns the full- 
grown cores of giant planets may ha ve been formed (|Cham bers 
2006). In the core accretion model (Polla cket al.lll996l) of gi- 
ant planet formation these cores of several M ffi attract significant 
amounts of gas from the nebula which results in the end in the 
birth of a gas giant planet. 

The alternative for the core accretion model is the gravi- 
tational instability scenario, in which a massive disk becomes 
gravitational unstable an d fragments into giant planets (iBossI 
1997; Maveretal. 2004). However, whether a disk can be- 
come unstable depends s ensitively on the cooling time scale 
(Pickett et al. 2000, 2003) and it is not yet clear if realistic pro- 
toplanetary disks will be subject to gravitational instabilities. 
Furthermore, the fragments produced in numerical simulations 
have masses around 10 Jupiter masses (Mj), which is rather high 
compared to the giant planets found in extrasolar planet surveys. 

One of the key differences between the two scenarios is the 
presence of a solid core. Unfortunately, the existence of solid 
cores in gas giant planets is still subj ect of debat e. Models of 
the internal structure of Jupiter ( Guillo t~t al.l2004l) can not con- 
strain the core mass from below due to uncertainties in the equa- 
tion of state near the center of the planet. However, the same 
models do suggest that Jupiter is enriched in solids with respect 
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to the Sun. If Jupiter was formed by means of gravitational in- 
stability it would have had to accrete these solids in a later stage. 

Accretion of solids also determines the luminosity of form- 
ing planets, which is important for future observations of pro- 
toplanetary disks. On their way to the core the accreted plan- 
etesimals deposit their energy in the gaseous envelope, which 
radiates away part of this energy into space. This accretion rate 
plays a major role in the dynamics of forming giant planets 
ilPollack etalll996l) . 

It is therefore important to un derstand the process of accre- 
tion of solids onto forming planets. Greenbe rg et alJ (11978) stud- 
ied accretion onto solid cores, and the effect of gas drag was in- 
vestigated by Weid enschilling & Davisl d!985h . who found that 
inward moving planetesimals will be captured in Mean Motion 
Resonances (MMRs). Only if the gas drag is strong enough (i.e. 
the planetesimals are small enough) they will make it to the sur- 
face of the planet. At first sight one might think that this slows 
down planetesimal accretion severely, but since their eccentric- 
ity is pumped up in the resonances collisions become much more 
frequent and the small debris resulting of catastrophi c collisions 
subsequently accrete onto the planet ( Weidenschilling & Davis 
119851) . However, iKarv et a"fl i 19931) showed that even particles 
small enough to move through all resonances do not always 
reach the surface of the planet. Instead, they are transferred into 
inferior orbits and they continue to move towards the central star. 

These studies including gas drag assumed a uniform gas 
disk, but hydrodynamical simulations of planets embedded in 
a gaseous disk show that planets more massive than 0.1 Mt 
start t o restructure their environment (e.g. D' Angelo et al. 2002, 
120031) . Eventually a deep gap forms in the disk for the most 
massi ve planets M p > Mj ( iLin & Papaloizoull993tlBrvden et all 
fl999l The pressure gradients associated with these gaps dramat- 
ically change the behavior of dust part icles, which may lead to 
a dus t gap while there is no gas gap ( Paardekooper & Mellema 
120041) . and also to a substantially redu ced accretion rate of 
solids dPaardekooper & Mellemal |2006al) . However, the latter 
study was limited in particle size and planet mass due to the 
fluid nature of the simulated dust component. In this paper we 
use a particle-based method to study dust accretion onto high- 
mass planets for arbitrary particle sizes. We focus on planets 
with masses large enough to affect their direct environment, 
M p > 0.1 Mj . This also allows us to do two-dimensional simu- 
lations, since the measured gas ac cretion rates for these planets 
in 2D are similar to the 3D values SD' Angelo et al.l2003l) . 

We start in Sect.FJJby reviewing the governing equations of 
motion, including the adopted gas drag law. In Sect. [3] we de- 
scribe the numerical method used to integrate the equations of 
motion, and Sect.0]is devoted to the disk model. In Sect. [5] we 
describe the results of the simulations, which we briefly discuss 
in Sect. [5] We conclude in Sect.Q 



2. Equations of motion 

Throughout we will work in a cylindrical coordinate frame (r, <p) 
with the central star in the origin. Because of the presence of a 
planet this is not an inertial frame, of which good use has been 
made in finding extrasolar planets, but we have to correct for this 
in the potential. The planet's orbit is circular, and the coordinate 
frame corotates with the planet at angular velocity Q p . Our unit 
of distance is the orbital radius of the planet, which is then lo- 
cated at (r, (p) = ( 1 , ti) throughout the simulation. The mass of the 
planet only enters the problem as a fraction of the stellar mass 
q = Mp/M*. When quoting explicit planetary masses (relative 



to the mass of Jupiter, Mj) we assume that the central star has a 
mass of 1 M . Then 1 Mj corresponds to q = 10~ 3 . 

2.1. Gas 

The evolution of the gas component of the disk is governed by 
the Euler equation s, which for this sp ecific case are described 
in detail in Paardekoop er & Mellemal J2006bl) . We do not solve 
the energy equation, but we use a locally isothermal equation of 
state: 



(1) 



where E g is the surface density of the gas, p is the vertically 
integrated pressure, and the isothermal sound speed c s is directly 
related to the disk thickness H: 



c s =HD. K = hvK 



(2) 



where Ok is the Keplerian angular velocity, h — H/r and vk = 
tQk- The gas has a kinematic viscosity v, which we take to be 
constant throughout the computational domain. We neglect the 
self-gravity of the gas. 



2.2. Dust Particles 

The equations of motion for a dust particle read: 

L 2 3<D 

r 3 dr 

<9<D 



dv 
dt 



+ fr 



dL 
dt 



(3) 
(4) 



where v is denotes radial velocity, L is the specific angular mo- 
mentum, <£ is the gravitational potential and / is the drag force, 
which we specify below. The potential contains contributions 
from the central star, the planet and indirect terms due to the 
acceleration of the coordinate frame. We do not consider self- 
gravity for the dust particles. 

2.3. Drag Force 

The nature of the friction between gas and dust depends on the 
size of the dust particles relative to the mean free path of the gas 
molecules. This is expressed by the Knudsen number: 



A 

Kn = — 

2s 



(5) 



where A is the mean free path of the gas molecules and s is the 
size of the dust particles under consideration. When Kn » 1 we 
are in the Epstein regi me of fr ee molecular flow, and the drag 
force is given by lSchaail963l) : 



F ep s = -7TS 2 p g |Av| Av 
1 1 



1 + 



4m A 



erf (to) + 



1 1 
— + 



2m 3 / yfn 



(6) 



where p g is the gas density, Av is the velocity of the dust particle 
relative to the gas, and to is the relative Mach number of the flow: 
to = |Av|/c s . Equation l|6} has the following asymptotic behavior: 



eps 



-^ps 2 PgCsAv 
-ns p g |Av|Av 



if |Av| <k c s ; 
if |Av| » c s . 



(7) 
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Because Eq. is difficult to implement numerically we inter- 
polate between thes e limits in or der to obtain the drag force for 
arbitrary velocities (Kwo kl 19751) : 



F - - 

* eps 



where 
/d = 



VT28^ 2 

— ~ — s p g c s f D \v 



r 9n 2 
1 H m 1 

128 



(8) 



(9) 



When Kn <K 1 the drag force is given by Stokes friction: 
F sto = -67r^ Dy u kin Av (10) 

where fi^m is the kinematic viscosity of the gas, which we write 

as: 



1 

/"kin = gPgVth'* 



(11) 



where v t h = V8 /nc s is the mean thermal velocity of the gas. Note 
that yUkin, and therefore also F st0 , is independent of the ambient 
gas density. This in contrast with F e ps, for which a higher gas 
density leads to a larger force of friction. 

In Eq. fllOi the drag coefficient ku depends on the Reynolds 
number Re as follows: 



( 1 +0.15/?e a687 
k D = j 3.96 1(T 6 Re 24 
O.URe 



if Re < 500; 

if 500 < Re < 1500; 

if Re > 1500. 



where 



Re 



2sp g \\v\ 

jUkin 



„ 7i m 



(12) 



(13) 



For flows of intermediate Knudsen number reliable expres- 
sions for the drag force are not readil y available, so we just in- 
terpolate (see lWoitke & Helling|2003l) : 



F = | 



' 3Kn 
. 3Kn + 1 



F f' 1 

eps \3Kn + l 



(14) 



Using Eqs. (JSJ, fllOi and (lilt we can rewrite the total drag force 

as: 



• V128tt 3K " ^ D + s 2 P r,c,KnAv 
(3Kn + l) 2 fe ' 



(15) 



Upon dividing by the mass of the dust particle, |^s 3 p p , where 
p p is the internal density of a dust particle, we can write for the 
force per unit mass: 



/ = pAv 

J T, 



(16) 



where Q. K is the Keplerian angular velocity and T s is the dimen- 
sionless stopping time: 



T, = J- 



(3Kn+ l) 2 



s Pp vk 

8 9Kn 2 f D + 3Kn k D r p g c s 



(17) 



Note that in the limit |Av| <k c s , Sc« 1 Eq. J 16b is linear in the 
relative velocity of the particle, which has great advantages for 
numerically solving the equations of motion (see Sect.0. 



3. Numerical method 

3.1. Gas 

We sol ve the flow equations for the ga s using the RODEO 
method ( Paardekooper & Mellema 2006b). This method was ex- 
tensively tested on the plan et-disk problem, and al so used in 
two-fluid mode in Paardekooper & Mellema (2004). It makes 
use of an approxi mate Riemann solver, togeth er with station- 
ary extrapolation (Eulderink& Mellema 1995) to integrate the 
source terms. See Paardekooper & Mellema (2006b) for details. 

3.2. Dust 

The dust component of astrophysical fluids can in principal be 
described as being a continuous fluid or as a collection of parti- 
cles. Both approaches have their advantages and disadvantages. 
The major advantage for the fluid approach is that the flow is 
accurately described even in regions of very low dust density. In 
the particle approach these regions would be severely depleted of 
particles, and therefore the local resolution would be relatively 
poor. In the cont ext of embed ded plan ets, the op ening of a dust 
gap in a gas disk ( Paardekooper & Mellema 2004 ) should prefer- 
ably be modeled using a dust fluid. 

How ever, the fl uid approach is not valid in all circumstances 
( iGaraud et alJ2 004). First of all, one fluid element (one grid cell) 
should contain enough dust particles to define a density, and to 
work with averaged velocities. Second, there should be enough 
interparticle collisions to make such an average meaningful. In 
a gas-dominated disk this means that the particles should couple 
reasonably well to the gas. In view of both these limitations of 
the fluid approach the larger particles of the dust size distribu- 
tion can not accurately be modeled as a fluid. The critical size 
depends on the gas density and the shape of the particles, but for 
dust grains larger than approximately 10 cm the fluid approach 
breaks down at 5 AU in the Minimum Mass Solar Nebula. If lo- 
cally the gas density is lowered severely, for example by a gap- 
opening planet, this maximum size goes down correspondingly. 

Therefore in cases of low gas density the particle approach 
seems to be the better option, although one has to make sure that 
enough particles reside inside the gap in order to obtain accurate 
results. In this study we focus on high-mass planets that open 
up gas gaps in the disk, and therefore we chose for the particle 
approach. 

In order to integrate the equations of motion numerically we 
used a standard second order symplectic integrator (leapfrog), in 
which alternatingly positions and velocities are updated. During 
a position update, velocities are assumed to be constant and vice 
versa: 



dr 

dt 
dv 

— = a(r,v), 
dt 



v constant 



r constant 



(18) 
(19) 



where a denotes the acceleration of the particle. 

The position update is trivial. For the radial velocity update 
we need to solve 



dv _ L 2 <9<D Qk a 



dt 



dr 



(20) 



Because we keep the gas velocity constant, the same differential 
equation applies to Av. In order to make the integration scheme 
suited in the regime r s « 1 we treat the drag term implicitly. 
Formally this can only be done in the limit fo, kr, = 1 (subsonic 
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analytical result with the numerical solution for T s — 0.005 in 
Fig-E and it turns out that they are indistinguishable. This shows 
that our implicit integration scheme for the drag force works cor- 
rectly. 

We model dust accretion onto the planet by taking away all 
particles within one half of the Roche lobe Rr of the planet: 



Rr 



(!)' 



(25) 



where q is the ratio of the masses of the planet and the cen- 
tral star. The actual solid core is much smaller than this, so we 
do not measure the number of impacts onto this core directly. 
However, the high-mass planets that we consider in this pa- 
per are able to capture a dense gaseous envelope that is able 
to efficien tly direct particles towards the center of the planet 
(see also Paardekooper & Mellema 2006a). We do not take away 
gas from the computational d omain, as has been done before to 
mimic gas accretion (iLubow et al.lll999t iD'Aneelo et alJ l2002 ; 
Paardekooper & Mellema 2006b). The process of accretion only 
influences the density structure close to the planet, and therefore 
this does not influence the results on dust accretion. 



Fig. 1. Radial velocity difference between gas and dust particles 
with T s = 0.005. 



and laminar regime). However, because supersonic drift veloci- 
ties as well as high Reynolds numbers only occur when r s » 1 
this is not a problem, because then the effects of gas drag are 
small and there is no need to treat the drag force implicitly. 
Upon differencing Eq. ( 1201 we obtain: 



Av - Av 



5t r> dr T s 



-Qk 



Av 



which has the solution 



Av 



1 + 



n K st 



(21) 



(22) 



A similar equation can be derived for the angular momentum 
equation. These equations are used to update the radial and an- 
gular velocities of the dust particles. 

The method was tested in two limiting cases: r s = oo and 
T s <K 1. When there is no coupling to the gas (T s = oo) the 
problem reduces to the restricted three-body problem, which has 
the Jacobi constant as an integral of motion: 



J = E - £l n ■ L 



(23) 



where E is the total energy of the particle under consideration. 
We reduced the magnitude of the time step for the particle inte- 
gration with respect to the hydrodynamical time step until / was 
conserved to 1 part in 10 6 after 100 orbits of the planet. This re- 
quired 10 particle time steps per hydrodynamical time step. We 
have adopted this time step for all simulations presented in this 
paper. 

In an axisym metric disk without a pl anet, the radial drift ve- 
locity is given bv ( iTakeuchi & Linl2002l) : 



Av = -t } (t; 1 +T s )v K 



(24) 



where r\ is the ratio of the gas pressure gradient to the stellar 
gravity. For uniform gas surface density, together with the tem- 
perature profile dictated by Eq. 0, r\ = h 2 . We compare this 



4. Initial and Boundary Conditions 

The ga s disk model is the same as in lPaardekooper & Mellemal 
(2004). The surface density of gas and dust particles is uniform 
initially, and we take the gas disk thickness to be H - 0.05 r, 
so h — 0.05, which is the canonical value for simulations of 
planet-disk interaction. In order to compensate for the radial 
pressure gradient dictated by Eq. Q the gas orbits at a slightly 
sub-Keplerian angular velocity, while the dust particles are on 
exact Keplerian orbits initially. We take all initial radial ve- 
locities of gas and dust to be zero. The gas has an anoma- 
lous viscosity v that is para metrized by the usual a-prescription 
dShakura & Sunvaevll973l) : 



a c s H 



(26) 



We adopt a constant value of v = 10~ 5 , which corresponds to 
a = 0.004 at the location of the planet. 

In order to calculate the stopping time in Eq. d!7i we need 
a value for the gas density. We fix the initial density to be 10" 11 
g crrT 3 at r = 1, which is appropriate for the location of Jupiter 
in the Minimum Mass Solar Nebula. From the gas density we 
can calculate the mean free path of gas molecules: 



mn 2 
np g r H2 



(27) 



where wh 2 and ru 2 are the mass and radius of an H2 molecule, re- 
spectively. The internal particle density p p is 1.25 g cm" 3 . Using 
these parameters, we can calculate T s as a function of particle 
size s (but note that r s also depends on |Av| through the factors 
fu and k D ). 

The dust density does not enter the equations of motion as 
long as it remains much smaller than the gas density. If not, the 
gas is affected by the drag force from the dust particles. We as- 
sume that initially the dust-to-gas ratio is equal to 1:100, and 
that we therefore can ignore the feedback on the gas. We vary 
the number of dust particles N typically from 5000 to 10 6 . 

The computational domain extends from r = 0.4 to r = 2.5 
in the radial direction, and from <p = to (f> — 2n in azimuth. This 
domain is covered by a uniform grid consisting of 128 radial and 
384 azimuthal cells which is used to evolve the gas component 
of the disk. 
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Fig. 2. Density structure after 50 orbits of a 1 Mj planet around a solar mass star. Left panel: particle distribution (T s = 0),N = 5000. 
Middle panel: particle distribution (T s =0),N— 10000. Right panel: logarithm of the gas surface density. 
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Fig. 3. Accretion of perfectly coupled particles onto a 1 Mj planet for four different amounts of particles 
mass in units of the total dust mass in the disk.Right panel: corresponding accretion rates, in disk masses 
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We ta ke the boundary c onditions for the gas to b e non- 
reflecting ( Godon 1996; Paardekooner & Mellema 2006b), in or- 
der to keep the waves generated by the planet from reflecting 
back into the computational domain. When a dust particle moves 
off the computational domain we remove it from the simulation. 



5. Results 

We have performed simulations with different particle sizes in 
order to sample the parameter space in stopping time ranging 
from the limit of perfect coupling T s = to the limit of no cou- 
pling T s - oo. The former limit offers an opportunity to compare 



the measured accretion rates to the gas accretion, so we will start 
by discussing perfectly coupled particles. 

When we increase the size of the particles the stopping time 
increases as well. We can distinguish three regimes: 

- r s « 1: the motion of the particles is dominated by gas drag. 
However, significant gas-dust separation may occur in the 
pres ence of planets on time scales o f approximately Q K 'r s 7 1 
(see Paardekooper & Mellema 2004). 

- T s » 1 : particle motion is dominated by gravitational inter- 
action with the star and the planet, with the gas drag as a 
small perturbation. 
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- T s » 1 : the orbital time scale for a particle is comparable to 
the time scale for gas drag. Depending on the gas density, 
there may be region in which either gas drag or gravity will 
dominate. 

5.1. Perfect coupling 

First of all we discuss the case for T s <K QkA?, where At is the 
magnitude of the time step as required for the integration of the 
Euler equations. For such a small value of the stopping time the 
dust particles are forced to move with the gas velocity, and the 
dust-to-gas ratio remains the same everywhere. In particular, gas 
and dust accretion rates should be equal as well, when we correct 
for the factor 100 in the dust-to-gas ratio. 

Gas accretion onto embedded planets h as been studied be- 
fore using different numerical approaches ([ Lutowej^alJ I 19991 
iD'Angelo etafl l2002t lPaardekooper~& Mellema 2006b), and 
they were found to agree reasonably well with each other. There 
is no strong dependence on the details of the accretion procedure 
(differences remain within ~ 40 %), except for intermediate mass 
planets for which the Roche lobe is approximately equal to the 
disk scale height. For these planets, accretion depends strongly 
on the conditions inside the Roche lobe and the measured accre- 
tion rates differ up to a factor of 2. Therefore we focus on a 1 
Mj planet, for which a well-defined gas accretion rate has been 
measured of M = 10~ 4 disk masses per orbit; a value that we 
would like to reproduce for the perfectly coupled dust particles. 
Typically this value was reached within 200 orbits of the planet. 

We have varied the number of particles in the simulation 
from 5000 to 10 5 . Note that in order to accurately measure an 
accretion rate the number of particles that is accreted onto the 
planet every orbit should be larger than one. For an accretion 
rate of 10~ 4 disk masses per orbit the total number of particles 
in the simulation N therefore needs to be larger than 10 4 . For a 
smaller value of N we need to rebin the accreted particles into 
larger time intervals. The loss of time resolution is of no impor- 
tance for the final accretion rate because it does not vary rapidly 
after ~ 50 orbits. 

In Fig. |2 we show the particle distribution for a simulation 
of a 1 Mj planet after 50 orbits. In the gas density (right panel) 
we see the gap that is developing, as well as the prominent spiral 
waves excited by the planet. For a run with 10000 particles (mid- 
dle panel) we can clearly see the gap and also the spiral waves in 
regions where the particle density is high enough. However, in 
the gap region near the planet, where in the gas density we see 
the origin of the spiral waves, nothing is to be seen in the particle 
density. This is because the spatial resolution in a particle simu- 
lation critically depends on the particle density. In regions of low 
density, as in the case of a disk gap, the resolution is correspond- 
ingly low and the spiral waves are not visible. This is of no con- 
cern for the outcome of these simulations, however, because the 
dust density has no dynamical effect. If the density was used for 
example in calculating the number of collisions between these 
particles, which in turn would create smaller dust grains, reso- 
lution effects would play an important role. The same would be 
true if the dust density becomes high enough that the gas feels 
its drag. The only thing we need to worry about for these sim- 
ulations is that we accrete enough particles per time step (see 
above). 

Even if we lower the number of particles N to 5000 there 
is still density structure visible (left panel of Fig. [2}- The outer 
edge of the gap is less sharp than in the case of N = 10000, but 
all features in the density can still be identified. 
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Fig. 4. Particle distribution (N = 10000) for T s = oo after 500 
orbits of a 1 Mj planet. Thick contours mark the boundary of 
the accretion region of the planet. The thin contour gives the 
corresponding zero-velocity curve. 



Also, if we look at the accretion rates for different values of 
N in Fig. [5] the results look promising, as expected. After 200 
orbits, the amount of disk mass that was accreted for the various 
number of simulated particles differs less than 15 %. Even more, 
the final accretion rates agree very well, yielding a mean value of 
M = 9.5 10" 5 disk masses per orbit, in close agreement with the 
value that came out of gasdynamical simulations. The total mass 
that is accreted as well as the slightly lower accretion rates can be 
attributed to the large imposed accretion inside the Roche lobe. 
While for the gas case one usually takes out a small fract ion of 
the density away from the accretion region every time step, for 
our particle simulations we remove all particles from the inner 
half of the Roche lobe. Such a high imposed accretion rate tends 
to decrease the total amount of ma ss that is accreted (see Fig. 7 
of Paardekooper & Mellema 2006b). 

We conclude that we are able to reproduce the gas accretion 
rates in the limit of perfectly coupled particles, and that we only 
need a relatively small amount of particles to achieve this. 

5.2. No coupling 

In this section we will focus on the limit of no coupling to the 
gas. For these particles the problem is equivalent to the restricted 
three-body problem, and therefore we can gain some insight by 
examining the integral of motion of the problem, J (see Eq. 
(I23» . In particular, the motion of a particle is confined to a re- 
gion bound by so-called zero-velocity curves, which can be ob- 
tained by solving Eq. i23\ for a given value of J — Jo and all 
velocities equal to zero. Particles with a starting value J = Jo 
will never cross the zero-velocity curve associated with /q. This 
means that the region in the disk from which the planet may 
accrete particles is also bounded, and once this feeding zone is 
empty accretion will stop. 
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Fig. 5. Accretion of uncoupled particles onto a 1 Mj planet for 
four different amounts of particles. 

In Fig. @] we show the particle distribution after 500 orbits 
of the planet for T s — oo and N = 10000. The region close to 
the planet is almost empty, and two annular empty rings appear 
at approximately r — 1.1 and r = 0.9. Particles in these rings 
travel on horseshoe orbits and suffer close encounters with the 
planet, which leads to accretion of these particles. Material that 
orbits further away from the planet (r > 1.1 and r < 0.9) can 
not reach the planet because of the constraint set by the Jacobi 
integral. The thick contours mark the radii at which particles may 
just make it to the planet, which is shown by the zero-velocity 
curve (thin contour) corresponding to the Jacobi constant at the 
location of the thick contour. 

The thick contours trace the outer edges of the feeding zone 
reasonably well, and almost all particles initially located just in- 
side this feeding zone are accreted by the planet. However, par- 
ticles closer to r — 1 do not suffer close encounters with the 
planet: they move on stable tadpole orbits and will not be ac- 
creted. In the solar system similar objects can be found around 
the orbit of Jupiter: the Trojan asteroids. The end result of the 
simulation is a particle distribution with two empty "rails" at a 
radial distance of 0. 1 from the planet, together with a non-empty 
corotation region. The empty semi-circle outside the orbit of the 
planet is due to resonant perturbations of the 2: 1 outer mean mo- 
tion resonance. 

In Fig. [5] we show the accretion rate up to 500 orbits for 
the same simulation, for four different values of N. Before t = 
100 (not shown in the figure) the accretion rate is very high, and 
the planet accretes roughly 10 % of the total disk during that 
time. However, after approximately 120 orbits the accretion rate 
drops below the gas accretion rate (10~ 4 disk masses per orbit) 
and it continues to decline for the whole simulation time. After 
500 orbits the particle accretion rate is an order of magnitude 
below the gas accretion rate, which makes further dust accretion 
negligible. 

All particles were started on circular orbits until now, while 
mutual gravitational inter actions between larger boulders lead s 
to eccentricity excitation dGreenzweig & LissauerlfT990l 119921) . 
The eccentricity of the particles' orbits determines the width of 



the accretion zone, and therefore the total mass that may be ac- 
creted. However, because the feeding zone is still bounded the 
accretion rate will still tend to zero after a certain amount of 
time. This is illustrated in Fig. [6] where we show the accre- 
tion rates for three different initial eccentricity distributions. All 
eccentricities in the specified range were uniformly distributed 
amongst the particles. It is clear that from approximately 130 
orbits the accretion rate is the same for all initial eccentricities, 
which is the same time scale as for the clearing of the feeding 
zone (see above). Therefore a non-zero initial eccentricity for 
the planetesimals enlarges the accretion zone, and therefore in- 
creases the total amount of mass that can be accreted (see the left 
panel of Fig. [6}, but it does not affect the time scale on which the 
feeding zone is cleared. 



- 5.3. Steady gas disk 



Clearly, the feeding zone needs to be replenished one way or the 
other in order for accretion to continue. In this section, we ex- 
plore the possibility that gas drag may prevent the feeding zone 
from getting empty. We do this by keeping the gas disk steady, 
in order to separate the effects on accretion of a non-uniform 
gas disk from effects due to pure dust dynamics. Als o this gives 
us an opportunity to compare with previous work bv lKarv et alJ 
(1993). 

Based on Eq. J24l > one might expect that dust accretion would 
always be stronger than gas accretion. Because the mass that 
flows inward through a circle with radius r equals: 



M = -In rlv r (28) 
the ratio of gas and dust accretion rates would be given by: 



M„ 



2tt rE g v r ,g 



= ,/|l + ^) 



(29) 



where d is the dust-to-gas ratio and Av is given by Eq. ( 1241 . In 
an unperturbed disk with constant surface density and kinematic 
viscosity v the radial velocity of the gas is equal to v r>g = — |r, 
and therefore we have: 



Av h 2 (T; 1 +T s )v K 2(r s -' + r s ) 



jac s h 



3 a 



(30) 



When all dust grains have the same stopping time T s the gas-dust 
mixture that is accreted by the planet has an effective dust-to-gas 
ratio of 



d* = d 



2{t; 1 +tX 

1 + — — 

3a 



(31) 



Using a = 0.004, already for T s = 0.001 the planet would get 
enriched in solids by more than 15 %. 

However, the planet m ay not be able to grab all the ma- 
terial that is passing by dKarv et all fl99l which leads to 
lower dust accretion rates. Gas, however, may be accreted at a 
higher rate than predicted by Eq. (12 8 1 ("see iLubow et alJ l 1999: 
tLubow & D'Angelo 200g). It is therefore not clear, even in an 
unperturbed gas disk, how dust accretion compares to gas accre- 
tion. 

In Fig. we show the dust accretion rates of dust particles 
onto a 1 Mj planet as a function of particle size. Note that we 
do not take into account effects of the particle size distribution 
in the disk: the accretion rates shown would be valid if all par- 
ticles in the disk would have a single size. The solid line gives 
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Fig. 6. Accretion of uncoupled particles (N - 5000) onto a 1 Mj planet for three different initial eccentricity distributions. Left 
panel: accreted dust mass in units of the total dust mass in the disk. Right panel: corresponding accretion rates, in disk masses per 
orbit. 



the accretion rate predicted by Eq. J28l > with the radial velocity 
of Eq. ( 1241 . Towards the lowest values of s (and therefore the 
lowest values of T s ) it approaches the value of the gas accretion. 
Note, however, that this is not the measured gas accretion rate, 
which is higher due to the high accretion efficiency of the planet 
dLubowetalll999l) . 

Particles r s « 1 have a small radial velocity, and they can 
essentially all be captured by the planet. This results in Fig. 
in accretion rates close to the predicted value for the smallest 
particles. When T s increases towards 1, the radial velocity of the 
dust particles increases, but at the same time is is easier for the 
planet to accrete them because they suffer no strong gas friction. 
Therefore essentially all particles with T s < 1 can be accreted by 
the planet. 

For larger particles a dramatic change in accretion rate oc- 
curs. Increasing the size of the particles by less than a factor 
of three leads to a reduction of accretion by almost two or- 
ders of magnitude. The reason for this low accretion rate lies in 
the m echanism of resonance trapping ( Wei denschilling & Davisl 
1985), where inward moving particles are captured in mean mo- 
tion resonances (MMRs). The drag force that moves the parti- 
cles inward is exactly balanced by resonant perturbations which 
are directed away from the planet. This is illustrated in Fig. |SJ 
where we show azimuthally averaged particle surface density as 
a function of particle semi-major axis a. Because the eccentrici- 
ties remain below 0.1 (see Fig. [5}, a almost coincides with r. 

The solid line in Fig. [8] corresponds to the particles size re- 
sponsible for the peak of the accretion rate in Fig.0 s = 37.0. 
There is a lot of structure near the position of the planet, but no 
gap is forming and there is no structure in the outer disk. Note 
that the surface density goes to zero at a « 2.25 because of the 
movement of the particles towards the central star. In the inner 
disk a small feature is visible at a — 0.48, which we can attribute 
to the 1 : 3 MMR. Note, however, that these structures in the 
inner disk are not stable: resonant perturbations as well as gas 
drag both push the particles inward. 
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Fig. 7. Accretion of particles onto a 1 Mj planet. The solid line 
gives the accretion rate of dust particles predicted by Eq. ( 1291 . 
All runs included 5000 particles. 



When we increase the particle size by less than a factor of 
three (dashed line in Fig. [8} the situation is remarkably different. 
Near a — 1 an annular gap starts to form, together with a strong 
overdensity at a — 1.2, which is close to the 4:3 mean motion 
resonance (a * 1.21). Apart from this large peak there is no 
structure in the outer disk. The inner disk shows a feature near 
the 1:2 MMR at a = 0.63. 



Sijme-Jan Paardekooper: Dust accretion onto high-mass planets 



9 



6 



4 



2 







- 




\ 










i 




T=1.0 










— — - T=2.7 


- 








- - T=27.0 


- 

!'i 




i 1 
i , 

1 

1 














; 


- ' 1 i i 




i ] 

w 


LA 














/-- 












\ / 


, , \.\ . ,\ 



0.5 1.0 1.5 2.0 

a 



Fig. 8. Particle surface density after 50 orbits of a 1 Mj planet 
and three different values for T s . 
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Fig. 9. Particle semi-major axis versus eccentricity after 200 or- 
bits of a 1 Mj planet and s = 1000 cm (T s * 27.0). 

Increasing the particle size by another factor of 10 again 
changes the picture (dash-dotted line in Fig. [3). Now the fea- 
ture at a = 1.2 has disappeared and instead we can see a peak 
near the 2:1 MMR at a — 1.6, together with a feature at a = 1.3, 
which corresponds to the 3:2 MMR. In the inner disk the 2:3 
MMR is excited at a = 0.76. 

It is clear from Fig.[8]fhat resonance trapping occurs when 
r s > 1, and that weaker-coupled particles may occupy lower- 
order resonances. Close to a resonance the eccentricity of a par- 
ticle is excited. In Fig.|9]we show the eccentricity e of the par- 
ticles' orbits for the T s = 27.0 case from Fig. [8] At the loca- 



tion of the 2:1 MMR e is excited to a maximum of approx- 
imately 0.05, which is consiste nt with the analytical result of 
Weidenschilling & Davis] ill 9851) who showed that the equilib- 
rium eccentricity for a particle in a (J + 1) : / MMR is approxi- 
mately given by e — 0.07 / -J j + 1 . 

Resonance trapping occurs because resonant perturbations 
tend to change the semi-major axis of the particle's orbit in such 
a way that it moves away from the planet. For outer resonances, 
this means that this effect opposes the drag-induced inward mi- 
gration. Inner resonances can never trap particles in a stable way, 
unless the drag force is directed outward. This may happen when 
a strong gas density gradient is present, possibly at the outer edge 
of a gas gap in the disk. 

The strength of a resonance (j + 1) : j depends on the mass 
planet-to-star mass ratio q as well as on the order of the reso- 
nance j. Higher order resonances, which are located closer to the 
planet, are stronger and they can stop inward moving particles 
more easily. Therefore the location at which a particle is stopped 
by resonant perturbations depends on the particle size s. This is 
clear from Fig. [S] where particles with T s = 27.0 (s = 1000 
cm) become trapped in the j — 1 resonance, while particles with 
r s = 2.7 (s - 100 cm) move all the way to the j — 3 resonance. 
Weidenschilling & Davis (1985) derived a minimum size for a 
particle that is able to move through the / h resonance: 

_ Pghr p 
3p p q C(j) f> 2 

where C(j) is an increasing function of j. For particles smaller 
than s m \ n the drag force is stronger than the resonant perturba- 
tions and the particles are able to move inward through the reso- 
nance. Because s min becomes smaller for higher values of j they 
may get trapped in a resonance closer to the planet. 

For higher values of j, successive resonances are more 
closely spaced and their overlap may lead to chaotic behavior 
(Wisdom 1980). The distance at which this happens is approxi- 
matelv dPuncan et alll989l) : 

\r-r p \ ~ 1.5 q 2p (33) 

For a planet of 1 Mj this means that particles cannot be trapped 
closer than \r - r p | = 0.2. It is interesting to note that this is 
also approximately equal to the width of the gas gap that is 
opened by a planet of the same mass (see iBrvden et alJll999l 
iPaarde kooper & Mellemal2006bl) . 

In Fig.[K)]we show s m i n as a function of j for a 1 Mj planet 
and a 0. 1 Mj planet. Also shown are the maximum values of j for 
which no chaotic orbits exist (from Eq. (I33». Focusing on the 1 
Mj planet we see that only particles smaller than approximately 
100 cm are able to travel through the j — 4 resonance, beyond 
which no stable trapping exists because of Eq. J33i . In terms of 
accretion onto the planet this means that only particles smaller 
than 100 cm are sufficiently coupled to the gas to make it all the 
way to the surface of the planet. From Fig.Qwe see that indeed 
the large jump in accretion rate happens around s = 100 cm. 
Our numerical result for s m ; n agrees within a factor of 2 with Eq. 

From Fig.Qit is clear that essentially all particles that make 
it through the resonances are eventually acc reted by the planet. 
This is not true in general (K arv et al] [l993). and especially for 
low-mass planets a significant fraction of particles that cross the 
orbit of the planet are transferred to the inner disk rather than 
to be accreted onto the planet. This is illustrated for the case 
of q = 10~ 4 in Fig. ^2 Moving from the smallest towards the 
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Fig. 10. Minimum particle size that can be trapped in the / h res- 
onance, for a 1 Mj planet and a 0.1 Mj planet. The vertical lines 
mark the onset of chaotic orbits for high j. 
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Fig. 11. Accretion of particles onto a 0. 1 Mj planet. The solid line 
largest particles (from left to right in Fig.[TTJ we see three classes gives the accretion rate of dust particles predicted by Eq. 
of particles: 



- I. Particles with s < 1 cm are very strongly coupled to the 
gas, which makes their radial velocity low enough for them 
to always reach the planet. 

- II. Particles with s > 1000 cm are again trapped in exterior 
resonances, and the accretion rate is therefore very low. 

- III. In between, there is a steady flow of particles onto the 
planet but there is a significant fraction of particles that 
misses the planet. 

The transition from Class III to Class II happens again approx- 
imately at the size predicted by Eq. J32i (see also Fig. I10i. The 
deviation from the accretion rate given by Eq. ( I29i for Class III 
particles seems to increase with s, except for the largest Class III 
particle at s = 100, which has a relatively high accretion rate. 

From Figs. an d EH we conclude that particles that are 
trapped in exterior resonances are excluded from accretion onto 
the planet. Only particles smaller than s m ; n may reach the planet, 
where s m j n is approximately given by Eq. J32I . Particles with for 
which 1 cm < s < s m [ n potentially have a higher accretion rate 
than the gas. 

5.4. Evolving gas disk 

However, when planets get massive enough to attract significant 
amounts of gas and dust they also start to reshape the surround- 
ing disk. Planets as massive as Jupiter open up a deep gas gap 
in the disk, but already Neptune-class planets restructure the gas 
disk in such a way th at a dust gap emerges, even tho ugh there 
is no gap in the gas ( Paardekoorjer & Mellema 2004). It was 
also show n that this leads to a drama tic decline in dust accre- 
tion ( Paardekoorjer & Mellema 2006a) for particles larger than 
0.1 cm. In this section we want to study the general case of dust 
accretion onto high-mass planets. 

We start by inspecting the particle distribution after 20 orbits 
of a 1 Mj planet in Fig.^] From the upper left panel to the lower 
right panel the stopping time increases from T s — (perfectly 
coupled particles, see Sect. 15. It to T s = oo (uncoupled particles, 



see Sect. 15. 21 . After 20 orbits the planet did not get a chance to 
open up a gap in the gas yet, as is clear from the upper left panel. 
The uncoupled particles in the lower right panel also do not clear 
a whole gap, but instead two empty "rails" (see also Fig.|4jl. 

For particles with T s = 0.027 (upper-middle panel in Fig. 
1 1 21 we see that the pressure gradients at the edges of the form- 
ing gap are pushing t he particles away from the planet's orbit, as 
was also observed in Paardekoorjer & Mellema (2004) for a ten 
times less massive planet. Because of the high mass of the planet 
in Fig. El an d the corresponding large pressure gradients, the 
effect takes place on much shorter time scales. Note also that be- 
cause of the forming gap the gas density is already significantly 
lower near the orbit of the planet, and consequently the particles 
couple not as well to the gas as at the start of the simulation. 
This, together with the larger pressure gradients at the gap edges 
accounts for the fast evolution of the dust component of 1 cm. 

When the stopping time approaches 1 the drift velocity of 
the particles becomes so large that already after 20 orbits of the 
planet an almost empty gap has formed (upper right and lower 
left panel of Fig. ll2b . The particles accumulate at the gap edges 
and at corotation. The latter because when a gas gap is form- 
ing the density at r = 1 is always larger than at |r — 1 1 = 0.1 
( see iKlevlfl 999t IPaa rdekoouer & Mellema 2006b), and the cor- 
responding pressure gradient pushes the dust particles towards 
r = 1. However, when the gas gap is sufficiently clean (which 
happens after approximately 100 orbits) this pressure gradient 
disappears and the particles near r = 1 will be dragged inward. 

When T s > 1 (lower middle panel in Fig. I12l > we see again 
that the pressure gradient associated with the edges of the gas 
gap stimulate the formation of the two empty "rails" compared 
to the case of no coupling. The first signs of resonant capture 
can be seen as a dark eccentric ring outside the planet's orbit. 
Because these particles do not couple well to the gas it becomes 
progressively harder to clean the corotation region while the gas 
gap is forming. Particles larger than 1000 cm are not removed 
from the corotation region before the gas gap is complete, which 



Sijme-Jan Paardekooper: Dust accretion onto high-mass planets 



11 




Fig. 12. Particle distribution after 20 orbits of a 1 Mj planet for 6 different particle sizes. Top row, from left to right: 5 = (perfect 
coupling, T s = 0), s = 1.0 cm (T s = 0.027) and s = 10.0 cm (T s = 0.27). Bottom row, from left to right: s = 37.0 cm (T s = 1.0), 
s = 1000.0 cm (T s = 27.0) and s — oo (no coupling, T s = oo). 



makes the horseshoe-shaped feature at r — 1 almost permanent 
just as in the case of no coupling. 

At later times, resonant trapping comes to dominate the den- 
sity structure of the largest particles. In Fig.^]we show the az- 
imuthally averaged surface density after 100 orbits of the planet 
for the gas and for three different particle sizes. The first thing 
that catches the eye is the large peak in the density of particles 
with s — 37 cm. The gas density (solid line in Fig. 1131 reveals 
that the pressure gradients around the peak all point away from 
r — 1 .4. This makes this location an efficient dust trap for basi- 
cally all material originally located at 1.1 < r < 1.8. Locally the 
surface density is increased with a factor of 25. 

The smaller particles (1 cm) move slower with respect to 
the gas than the 37 cm particles. As a result two peaks can be 
seen near r — 1 .4: one from outward moving particles coming 
from r « 1 and one from inward moving particles coming from 
r > 1.7. The large particles of 1000 cm also show a double- 
peaked feature. Note that the density gradient near the 2: 1 MMR 
at r = 1.6 is such that the particles are dragged inward more 
efficiently than in the case of the steady gas disk. Therefore while 
Fig.^Spredicts that the 1000 cm particles should be captured in 
the 2: 1 MMR they are able to move all the way to r ~ 1 .4. 

In Fig.[2]we show the time evolution of the surface density 
of particles of 37 cm. Four distinct features can be distinguished 
at r ~ 0.6, r ~ 1, r » 1.4 and r » 1.9, which are all due to 
the density structure in the disk. We have already mentioned the 
large peak near r — 1 .4, which builds up largely from the inside 



within 100 orbits. Note that the peak slowly moves outward but 
stops near r = 1 .4, because it has arrived at a pressure maximum. 

For the peak near the inner edge of the gas gap at r = 0.6 the 
situation is slightly different. Again, the peak is moving slowly 
but this time the density gradient inside r = 0.6 is not strong 
enough to stop the particles completely. The peak slows down, 
but continues steadily towards the inner boundary and in the end 
it will move off the grid. 

The peak at corotation is steady. Because the gas density has 
dropped approximately 2 orders of magnitude, the effective stop- 
ping time at r — 1 has increased from 1.0 to 100.0. This means 
that it will take several hundreds of orbits more to remove this 
feature. Finally, a small outward moving feature can be seen as- 
sociated with the density gradient near r — 1.9. Note that reso- 
nant trapping is not important for particles of s — 37 cm, which 
agrees with Fig. [7] However, the density gradients induced by 
the planet are so large that even particles of 1000 cm are not 
captured in the 2: 1 MMR. 

The behavior of particles of 1 cm is shown in Fig.^] First 
thing to note is the absence of the density feature near r = 1 .9, 
which is due to the stronger gas drag that makes the particles 
follow the inward viscous motion of the gas. Near r = 1.4 we 
can clearly see the two density peaks originating from both sides 
of the gas density maximum slowly moving towards each other. 
In contrast with the 37 cm particles from Fig.[3]the outer peak 
is the strongest, which is due to the lacking feature near r = 1.9. 
Particles of 37 cm originally located at r = 1.8 move outward, 



12 



25 



2(1 



10 



Sijme-Jan Paardekooper: Dust accretion onto high-mass planets 

10 I — i — 1 — 1 — 1 — 1 — i — 



gas, t=200 
s=l cm 
s=37 cm 
s= 1000 cm 




gas, t=200 

1=20 

1=50 

t=100 

t=200 




Fig. 13. Azimuthally averaged surface density of gas and dust Fig. 15. Particle surface density (s = 1 cm) around a 1 Mj planet 
after 100 orbits of a 1 Mj planet for three different particle sizes, at different times. 
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Fig. 16. Accretion of particles onto a 1 Mj planet and an evolv- 
ing gas disk for four different particle sizes. After 30 orbits, all 
Fig. 14. Particle surface density (s = 37 cm) around a 1 Mj accretion rates stay exactly zero, 
planet at different times. 



while their 1 cm counterparts move inward and contribute to the 
inward moving peak near r = 1 .4. 

From Fig. [2] we conclude that the particles that had the 
largest accretion rate in Fig. actually create the cleanest dust 
gaps. Already after 20 orbits the feeding zone is almost empty, 
and this has severe implications for the accretion rates. In fact, 
for all particles that open dust gaps the accretion rate equals zero 
after a finite amount of time. 

This is illustrated in Fig. where we show the accretion 
rate for four different particle sizes. After only 30 orbits, all ac- 



cretion rates stay exactly zero. Note that this is unlike all previ- 
ous models: for perfectly coupled particles as well as uncoupled 
particles there was still residual accretion after 200 orbits (see 
Fig- 13- This shows how powerful pressure structure in the disk 
is when it comes to moving dust particles. 

The time scale at which accretion stops depends on the stop- 
ping time. The shortest time scales correspond to particles of 
10 and 100 cm, which have an initial stopping time of T s * 1. 
Smaller particles as well as larger particles have a smaller drift 
velocity (see Eq. (I24» and therefore it takes longer for them to 
clear the dust gap. Eventually, for the smallest particle sizes we 
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Fig. 17. Accretion of particles onto a 1 Mj planet and an evolving Fig. 18. Accretion of particles onto a 0.1 Mj planet and an evolv- 
gas disk after 200 orbits. The solid line gives the gas accretion ing gas disk after 200 orbits. The solid line gives the gas accre- 



rate. All runs included 5000 particles. 



tion rate. All runs included 5000 particles. 



expect to recover the gas accretion rate, while for the largest par- 
ticle sizes we expect some residual accretion. 

The accretion rate during the first 5 orbits of the planet 
grows monotonically with particle size. Particles of 1 cm ac- 
crete approximately as fast as perfectly coupled particles in this 
stage (see Fig- EJ. Our resolution close to the planet is not high 
en ough to see the Roche lobe cleari ng stage as was observed 
in Paardekooper & Mellema (2006a) using an adaptive mesh. 
Larger particles accrete faster, and the particles of 1000 cm have 
an even higher accretion rate than uncoupled particles (see Fig. 
[6jl. The turnover occurs around T s = 1, which marks the tran- 
sition from gas drag dominated dynamics (T s < 1) and grav- 
ity dominated dynamics (T s > 1). During this stage, when the 
planet has not opened up a gas gap yet, gas drag is able to speed 
up accretion as long as we are in the gravity dominated regime. 
However, when the initial stopping time T s < 1 gas drag really 
takes over and forces the particles to accrete at the same rate as 
the gas. 

In Fig. we show the accretion rate after 200 orbits of a 
1 Mj planet for various particle sizes. Again we see that the ac- 
cretion rate is identically zero for a large range of sizes, while 
at both ends of the figure we approach the accretion rates of the 
limiting cases of perfectly coupled particles (left) and uncoupled 
particles (right). The dust accretion rates almost never become 
steady, unless they are identically zero or in the limiting case of 
perfectly coupled particles (cf. Figs. [3] and The radial distri- 
bution of the smallest particles evolves on a time scale propor- 
tional to 7"" 1 , which becomes prohibitively long when r s « 1. 
Therefore we chose a conservative approach to take the accre- 
tion rate after 200 orbits, when the gas accretion has settled to a 
steady state, as the final accretion rate for a given particle size, 
even though the particle accretion rate is still declining. 

The transition from no accretion to gas accretion begins 
around s = 0.001 cm, or equivalently T s = 0.000027. This 
means that the total accretion flow onto the planet is severely 
depleted in particles larger than 10 /itn. Even more, dust accre- 



tion never exceeds gas accretion, which makes it very difficult to 
enrich a giant planet in solids. 

A ten times less massive planet of 0. 1 Mj does not open up 
a gas gap, and therefore the planet-i nduced pressure grad ients 
in the d isk are less strong. However, Paardekooper & Mellema 
(2006a) showed that even for such a small planet dust accretion 
severely slows down. Figure ^] shows the dust accretion rates 
for various particle sizes. Again, we see accretion rates that are 
identically zero for s > 1 .0 cm, but the transition from gas accre- 
tion rates to zero accretion is shifted to larger s with respect to 
the planet of 1 Mj. We find that particles larger than 100 fim have 
significantly low dust accretio n rates compared to the gas, which 
agrees with the result of Paardekooper & Mellema ( 2006a). Also 
for this low-mass planet it is very hard to get enriched in solids. 

More massive planets than Jupiter open up wider and deeper 
gaps, but the gap edges do not get much steeper. Therefore the 
minimum particle size that is able to move across the outer gap 
edge to be accreted by the planet is the same as for the 1 Mj 
planet. The width gas gap as well as the dust gap increases, 
though, which brings the density peak at the outer gap edge close 
to the 2:1 MMR. This could lead to interaction between large 
bodies trapped in the 2:1 MMR and smaller bodies trapped at 
the gap edge. 

6. Discussion 

When IWeidenschilling & Davisl (Q985) discovered that orbital 
resonances are effective dust traps, they argued that this would 
not stop dust accretion completely. The increase in eccentric- 
ity near the resonance would promote collisions between plan- 
etesimals, and the resulting debris would be small enough to 
move through the resonances towards the planet. iKarv et all 
( 1993) showed that this is not the whole story: particles small 
enough to make it through the orbital resonances have a sig- 
nificant chance of being transferred to interior orbits rather 
than being accreted onto the planet. We have demonstrated that 
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there is yet another barrier for smaller dust particles when the 
planet is massive enough to induce significant pressure gradi- 
ents in the disk, which is already the case for a planet of 1 5 M ffl 
(Paardekooper & Mellema 2006a). 

According to Fig. [TO] when a planetesimal breaks up into 
100 cm sized boulders in a catastrophic collision, these can sub- 
sequently be accreted by a 1 Mj planet. However, inspection of 
Fig. ^[reveals that 100 cm boulders are stopped at the outer gap 
edge, and they need to be grinded down to 10-100 /mi before 
they can reach the planet. Moreover, it becomes progressively 
harder to break smaller particles because they are better coupled 
to the gas and their relative velocities are low. This means that 
the accretion flow onto the planet will be severely depleted in 
solids. The same holds for a 0.1 Mj planet, although somewhat 
larger particles may enter the accretion flow (see Fig. ll8l l. 

This has severe conse quences for t h e fina l composition of 
the planet. According to Pollack et alJ dl996l) . a 1 Mj planet 
at 5.2 AU approximately spends the final 500,000 years of its 
formation accreting the major part of its massive gaseous en- 
velope. Before that happened, the planet had already reached 
a mass of 0.1 Mj , which makes it vulnerable to dust gap for- 
mation. Therefore the major part of the mass accreted by the 
planet is relatively dust-poor compared to the interstellar value, 
which is of importance for the final composition of the planet. 
When Jupiter starts accreting gas when its solid core as well 
as the gaseous envelope both have a mass of 10 M ffi the differ- 
ence in solid content between dust-rich gas accretion and dust- 
poor gas accretion can be as large as 33 %. Therefore, the en- 
richment in solids in Jupiter requires traditional enrichment sce- 
narios (e. g. core erosion dGuillot et alJl2004l) . clathrate hydrate 
trapping (Gaut ier et al.l20 01a b)) should be much more effective 
than assumed until no w. 

d Patterson (1987), in cont rast with the suggestion by 
denschilling & Davis (1985), suggested that the accumula- 
tion of bodies in orbital resonances may lead to enhanced growth 
of planetesimals. Especially for a planet on an eccentric orbit, 
for which the eccentricities of the planetesimals in low-order 
resonances stay rather low. In either case, the size distribution 
of trapped particles may be very dynamical, because growth as 
well as destruction is enhanced with respect to the rest of the 
nebula. This makes it hard to predict the outcome of resonance 
trapping with respect to the resulting size distribution. 

Another interesting question is what the change in opacity 
does to the accretion flow. Dust is the major source of continuum 
radiation in the disk, and if the dust content is very low the gas 
may not be able to cool efficiently through radiation. This may 
severely slow down gas accretion. The magnitude of the opacity 
drop depends on how many of the smallest dust particles enter 
the accretion flow, which in turn depends on the details of the 
collisions of the larger bodies. 

We have focused on planets of relatively high mass, in order 
for the two- dimensional approach to be valid. For planets more 
massive than 0.1 Mj the gas accretion rates i n two and three 
dimensional simulations agree reasonably well ( iD'Angelo et al.l 
120031) . and because the dust particles are probably confined to 
an even thinner layer we believe that three-dimensional effects 
will be of minor importance. Also the results presented here do 
not depend on the density structure within the Roche lobe of the 
planet. 

Tu rbulence in ci rcumstellar disks is probably of magnetic 
origin (Balb us & HawievTl 990). If strong enough, it may prevent 
the formation of a deep gap but for reasonable disk parameters 
the g ap formation paradigm is not changed ( Papal oizou et al.l 
2004). Note that the only thing that is needed to trigger dust 



gap formation is the existence of planet-induced pressure gradi- 
ents, and therefore as long as the turbulent disk structure does 
not overcome gap formation the mechanism presented here will 
operate to reduce dust accretion onto high-mass planets. 

The effect of a nonzero eccentricit y of the p l anet's orbit 
on resonance trapping was studied by IPattersonI ( 119871) and 
iKarv & Lissauerl (JT995). The effect on the gas disk, in partic- 
ular the form ation of a gap, is not different from case of a 
circular orbit ( Papaloizou et al. 2001). Moreover, eccentricity 
growth of e mbedded pla nets is associated with deep annular gaps 
(lArtvmowiczlll992t [Lin & Papaloizoulll993l [Goldreich & Saril 
2003: ISari & Goldreichl2004UKlev & Dirksenl2006l) . Therefore 
it is likely that dust accretion onto eccentric planets also suffers 
from the planet-induced pressure gradients in the disk. 

In this paper, we have not considered migrating planets. Gap- 
opening planets experience Type II migration, for which the mi- 
gration speed is comparable to the viscous inward drift of the 
gas. This migration rate is not enough to change the results de- 
scribed in the previous section. However, for planets that open 
up a shallow gas gap, the inward migration speed can be much 
higher. It is then possible that, depending on their size, particles 
will approach the planet both from the inner disk and the outer 
disk. The planet will overtake the smallest particles, while the 
particles that are less well-coupled to the gas still overtake the 
planet on their way in. This means that dust trapping in the inner 
resonances becomes possible for particles that are being over- 
taken by the planet. Furthermore, the inner gap edge acts as a 
dust barrier for particles coming from the inner disk in the same 
way as the outer gap edge does for particles coming from the 
outer disk. Therefore, also in the case of a migrating planet dust 
accretion can be slowed down significantly. 

When a planet is slowly accreting gas and dust from the sur- 
rounding disk, it will very gradually approach the mass that is 
required to open up a gap, first only in the dust disk, later also in 
the gas disk. Before that time, only particles that are captured in 
resonances will not reach the planet. This applies to the largest 
particles. As the planet gains mass, it may capture smaller and 
smaller particles in resonances. Therefore, the maximum particle 
size that may accrete onto the planet is decreasing with planetary 
mass (see also Figs.ll7landll8>. 

For simplicity, we have chosen a constant initial surface den- 
sity. When the density profile is steeper, the surface density in 
the outer disk will be lower than in the case considered here. 
Therefore, particles in the outer disk will be less well coupled, 
the effects described in the previous section will be even more 
pronounced. 



7. Summary and conclusion 

In this paper we have studied the accretion rates of solid particles 
of different sizes onto high-mass planets using a new particle- 
based method for the integration of the equations of motion of 
the dust component. We have shown that we can reproduce the 
gas accretion rates as found in hydrodynamical simulations in 
the limit of perfectly coupled particles. In the limit of no cou- 
pling to the gas, the planet is able to clear its feeding zone in 
approximately 100 orbits, regardless of the initial eccentricity 
distribution of the particles. In order for accretion to proceed it 
is necessary to continuously replenish the feeding zone. 

For a non-evolving gas disk three different classes of parti- 
cles can be distinguished: the smallest particles couple very well 
to the gas and they accrete at the same rate as the gas, the largest 
particles get trapped in exterior resonances with the planet, and 
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only particles with intermediate sizes (1 cm < s < s m i n ) accrete Youdin, A. N. & Shu, F. H. 2002, ApJ, 580, 494 
onto the planet at a higher rate than the gas. 

When the gas is allowed to evolve dynamically it turns out 
that the planet-induced pressure gradients due to the formation 
of an annular surface density depression play a major role in the 
movement of dust particles towards the planet. Only the particles 
that couple very well to the gas are able to move across the outer 
edge of the density depression, which limits the size of particles 
that can be accreted by a 1 Mj planet to s < 10 /im. For planets 
of lower mass somewhat larger particles may be accreted, but 
even for a 0.1 Mj planet s < 100 //m. The lack of accreting large 
bodies may have important consequences for growing planets in 
disks, especially with regard to their enrichment in solids. If a 
large mass fraction of the solid component of the disk resides in 
particles larger than 10 yum then it is not possible to enrich the 
planet. 
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